##########################################
## "Voters get what they want"          ##
##########################################
## Heinrich, Kobayashi, Long            ##
##########################################

##          
## Auxiliary functions
## June 29, 2017


## Function that estimates interaction models
#############################################
estimate_A1_interaction <- function(to_estimate, ABCs, ...)
{
  this_dv <- "cbind(lnaidpc_l, lnaidpc) ~ 1 +"
  this_dv_tmp <- "lnaidpc ~ (1|countrynumcode_g) +"
  lthis_dv <- " + llnaidpc"
  
  this_form <- "polity2 + lnworldaidtotal + lln_rgdpc + lln_population +
  + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar +  
  + ldonorallyneighbor2 + ls3un + llnreftotal"
  
  if(A[to_estimate] == "Security hierarchy")
  {
    this_data <- data_SH
    this_form <- paste0("lphysint + us_SH1995 + llnNYT_leader + I(us_SH1995 * lphysint) + I(llnNYT_leader * lphysint) + 
                        I(llnNYT_leader * us_SH1995) + I(llnNYT_leader * us_SH1995 * lphysint) ", lthis_dv, " + ", this_form)
  }
  if(A[to_estimate] == "Economic hierarchy")
  {
    this_data <- data_EH
    this_form <- paste0("lphysint + us_EH1995 + llnNYT_leader + I(us_EH1995 * lphysint) + I(llnNYT_leader * lphysint) + 
                        I(llnNYT_leader * us_EH1995) + I(llnNYT_leader * us_EH1995 * lphysint) ", lthis_dv, " + ", this_form)
  }
  
  
  ## Estimate model
  #################
  out <- vector("list", chains)
  for(k in 1:chains)
  {  
    tmp <- MCMCglmm(fixed=as.formula(paste0(this_dv, this_form)),
                    random= ~ us(1):countrynumcode_g,
                    data=this_data, verbose=round(niter/5),
                    prior=list(R=list(V=1, nu=.001), G=list(G1=list(V=1, nu=.001))),
                    family="cengaussian", nitt=niter, burnin=burnin, thin=thin)
    out[[k]] <- tmp
  }
  out_tmp <- lmer(formula=as.formula(paste0(this_dv_tmp, this_form)),
                  data=this_data)
  
  
  ## Calculate marginal effects | hierarchy
  #########################################
  hier_grid <- seq(from=0, to=1, by=.1)
  beta <- ldply(.data=out, .fun=function(x) x$Sol)
  beta_here <- beta[, c(2, 5, 6, 8)] 
  news_tmp <- this_data$llnNYT_leader
  
  me_effs <- c()
  
  for(j in 1:length(hier_grid))
  {
    for(i in 1:nrow(news_sets))
    {
      ## Trim to proper news subset    
      news_here <- news_tmp[news_tmp >= quantile(news_tmp, p=news_sets$Lower[i]) & news_tmp <= quantile(news_tmp, p=news_sets$Upper[i])]
      
      ## Marginal effect of human rights violations
      me_physint <- as.matrix(beta_here) %*% t(cbind(1, hier_grid[j], news_here, hier_grid[j] * news_here))
      ## Marginalize over news
      me_physint <- rowMeans(me_physint) 
      q_s <- quantile(me_physint, probs=c(.025, .5, .975))       
      
      ## Marginal effect of marginal effect
      me_me_physint <- as.matrix(beta[, c(5, 8)] ) %*% t(cbind(1, news_here)) 
      ## Marginalize over news
      me_me_physint <- rowMeans(me_me_physint)
      
      me_effs <- rbind(me_effs,
                       data.frame(Hier=hier_grid[j],
                                  HierType=A[to_estimate],
                                  NewsSet=news_sets$Name[i], 
                                  NewsType="Leader",
                                  Q025=q_s[1], Q500=q_s[2], Q975=q_s[3],
                                  ProbNegative=mean(me_physint < 0),
                                  ProbMEChange=mean(me_me_physint > 0)))
    }
  }
  
  ## Convergence
  ##############
  post_list <- llply(.data=out, .fun=function(x) mcmc(cbind(x$Sol, x$VCV)))
  rhat <- max(gelman.diag(mcmc.list(post_list))$psrf[,1], decreasing=TRUE)
  load("output/Rhat.Rdata")
  Rhat <- rbind(Rhat, data.frame(Rhat=round(rhat, di=2), Which=paste0("Model A1-", to_estimate)))
  save(Rhat, file="output/Rhat.Rdata")
  
  
  ## Save output
  ##############
  output <- list(SubstEff=me_effs,
                 Model=out,
                 convergence=rhat,
                 lmer_model=out_tmp)      
  save(output, file=paste0("output/analysis1_interaction/Output-", to_estimate, ".Rdata"))
  
  return(rnorm(1))
}


## Function that estimates subset models
########################################
estimate_A1_subset <- function(to_estimate, ABCs, ...)
{
  this_dv <- "cbind(lnaidpc_l, lnaidpc) ~ 1 +"
  this_dv_tmp <- "lnaidpc ~ (1|countrynumcode_g) +"
  lthis_dv <- " + llnaidpc"
  
  this_form <- "polity2 + lnworldaidtotal + lln_rgdpc + lln_population +
  + lln_trade + lalliance + socialist + ColdWar + I(ColdWar * socialist) + lwar +  
  + ldonorallyneighbor2 + ls3un + llnNYT_leader + llnreftotal"
  
  ## Modification based on hierarchy measure
  if(ABCs$A[ABCs$FM$A[to_estimate]] == "Security hierarchy")
  {
    this_data <- data_SH
    this_form <- paste0("lphysint + us_SH1995 + I(us_SH1995 * lphysint) ", lthis_dv, " + ", this_form)
  }
  if(ABCs$A[ABCs$FM$A[to_estimate]] == "Economic hierarchy")
  {
    this_data <- data_EH
    this_form <- paste0("lphysint + us_EH1995 + I(us_EH1995 * lphysint) ", lthis_dv, " + ", this_form)
  }
  
  ## News subsetting
  if(ABCs$C[ABCs$FM$C[to_estimate]] == "Leader")
  {
    this_c_lower <- quantile(this_data$llnNYT_leader, ABCs$B$Lower[ABCs$FM$B[to_estimate]])
    this_c_upper <- quantile(this_data$llnNYT_leader, ABCs$B$Upper[ABCs$FM$B[to_estimate]])
    this_data <- subset(this_data, llnNYT_leader >= this_c_lower & llnNYT_leader <= this_c_upper)
    if(sd(this_data$llnNYT_leader) == 0) this_form <- str_replace(string=this_form, pattern="\\+ llnNYT_leader",
                                                                  replacement="")
  }
  if(ABCs$C[ABCs$FM$C[to_estimate]] == "Human rights")
  {
    this_form <- str_replace(string=this_form, pattern="llnNYT_leader", replacement="llnNYT_HR")
    this_c_lower <- quantile(this_data$llnNYT_HR, ABCs$B$Lower[ABCs$FM$B[to_estimate]])
    this_c_upper <- quantile(this_data$llnNYT_HR, ABCs$B$Upper[ABCs$FM$B[to_estimate]])
    this_data <- subset(this_data, llnNYT_HR >= this_c_lower & llnNYT_HR <= this_c_upper)
    if(sd(this_data$llnNYT_HR) == 0) this_form <- str_replace(string=this_form, pattern="\\+ llnNYT_HR",
                                                              replacement="")
  }
  
  ## Estimate model
  #################
  out <- vector("list", chains)
  for(k in 1:chains)
  {  
    tmp <- MCMCglmm(fixed=as.formula(paste0(this_dv, this_form)),
                    random= ~ us(1):countrynumcode_g,
                    data=this_data, verbose=round(niter/5),
                    prior=list(R=list(V=1, nu=.001), G=list(G1=list(V=1, nu=.001))),
                    family="cengaussian", nitt=niter, burnin=burnin, thin=thin)
    out[[k]] <- tmp
  }
  out_tmp <- lmer(formula=as.formula(paste0(this_dv_tmp, this_form)),
                  data=this_data)
  
  
  ## Calculate marginal effects | hierarchy
  #########################################
  hier_grid <- seq(from=0, to=1, by=.1)
  beta <- ldply(.data=out, .fun=function(x) x$Sol)
  
  me_physint <-  beta[, 2]  + beta[, 4] %*% t(hier_grid)
  q_s <- colQuantiles(me_physint, probs=c(.025, .5, .975))
  
  me_effs <- data.frame(Hier=hier_grid,
                        HierType=ABCs$A[ABCs$FM$A[to_estimate]],
                        NewsSet=ABCs$B$Name[ABCs$FM$B[to_estimate]], 
                        NewsType=ABCs$C[ABCs$FM$C[to_estimate]], 
                        Q025=q_s[,1], Q500=q_s[,2], Q975=q_s[,3], 
                        ProbNegative=colMeans(me_physint < 0),
                        ProbMEChange=mean(100 * beta[,4] > 0))
  
  
  ## Convergence
  ##############
  post_list <- llply(.data=out, .fun=function(x) mcmc(cbind(x$Sol, x$VCV)))
  rhat <- max(gelman.diag(mcmc.list(post_list))$psrf[,1], decreasing=TRUE)
  load("output/Rhat.Rdata")
  Rhat <- rbind(Rhat, data.frame(Rhat=round(rhat, di=2), Which=paste0("Model A1-", to_estimate)))
  save(Rhat, file="output/Rhat.Rdata")
  
  
  ## Save output
  ##############
  output <- list(SubstEff=me_effs,
                 Model=out,
                 convergence=rhat,
                 lmer_model=out_tmp)      
  save(output, file=paste0("output/analysis1_subset/Output-", to_estimate, ".Rdata"))
  
  return(rnorm(1))
}



## Log-posterior density of randomized response with
## four levels of the outcome
####################################################
rr.logpostdens.k4 <- function(theta, ...)
{
  beta <- theta[3:length(theta)]
  alpha <- c(1.2^(theta[1]),
             1.2^(theta[1]) + 1.2^(theta[2]))
  
  lp <- X %*% beta
  prZ1 <- pnorm(0 - lp)
  prZ2 <- pnorm(alpha[1] - lp) - pnorm(0 - lp)
  prZ3 <- pnorm(alpha[2] - lp) - pnorm(alpha[1] - lp)
  prZ4 <- 1 - pnorm(alpha[2] - lp)
  
  ll1 <- log(0.5 * prZ1 + 0.5 * 0.25)
  ll2 <- log(0.5 * prZ2 + 0.5 * 0.25) 
  ll3 <- log(0.5 * prZ3 + 0.5 * 0.25) 
  ll4 <- log(0.5 * prZ4 + 0.5 * 0.25) 
  
  ll <- sum(I(Y == 1) * ll1 + I(Y == 2) * ll2 + I(Y == 3) * ll3 + I(Y == 4) * ll4)
  pd <- ll + sum(log(dnorm(theta, 0, 5)))
  return(pd)
}


# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
# Source: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
###################################################################################
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


## Auxiliary function
#####################
substrRight <- function(x, n){
  substr(x, nchar(x)-n+1, nchar(x))
}



